Energy drift in reversible time integration

Energy drift is commonly observed in reversible integrations of systems of molecular dynamics. We show that this drift can be modelled as a diffusion and that the typical energy error after time T is O( √ T ). PACS numbers: 45.10.−b, 05.45.Pq In simulations of conservative systems, the energy H is usually monitored as a check on the calculation. In symplectic integration of Hamiltonian systems, it is known that the integrator is very close to the flow of a Hamiltonian system with Hamiltonian close to H, so that one can give conditions under which the energy error is bounded for exponentially long times [14]. However, in reversible integration [1–3, 5, 7, 8, 10], one typically sees the energy drift away from its initial value. In this letter we model this drift as a diffusion process, showing that the expected drift after time T is O( √ T ). There are several reasons why one might use a reversible integrator on a conservative system. First, if the system is Hamiltonian, a symplectic integrator might be prohibitively expensive. This occurs if one wants to adaptively vary the time step, which can be much cheaper to do reversibly than symplectically, or if the symplectic structure is noncanonical, perhaps as a result of a change of variables. See, e.g., the discussion of the Nosé–Hoover thermostat of molecular dynamics in [2]. Second, if the system is not Hamiltonian but still has a first integral H, then a reversible integrator is the natural choice of geometric method. One can construct integrators which are reversible and preserve energy, but they are expensive, typically fully implicit in the dependent variables and in the (introduced) Lagrange multipliers. It is usually much cheaper to preserve just the reversibility, which is the dominant property characterizing the dynamics, and merely monitor the energy. We consider systems with phase space M and dynamics ẋ = f (x), reversible under the diffeomorphism R : M → M , i.e. R∗f = −f or f (R(x)) = −TxR · f (x) ∀x ∈ M , and 0305-4470/04/450593+06$30.00 © 2004 IOP Publishing Ltd Printed in the UK L593 L594 Letter to the Editor with a symmetric first integral H : M → R, i.e. Ḣ = f (H) = 0 and H(R(x)) = H(x). (H can be any symmetric integral; we are calling it ‘energy’ for illustration.) The integrator, a diffeomorphism φτ : M → M , which approximates the time-τ flow of f , is also assumed to be reversible, i.e. φτ (R(x)) = R ( φ−1 τ (x) ) . We are interested in the energy drift: for a given initial condition x0 ∈ M , which we take for convenience to have zero energy, H(x0) = 0, what is the behaviour of the sequence { H ( φ τ (x0) )} ? Extensive numerical evidence suggests that it is not bounded but wanders erratically: in fact, it looks like a random walk. This is in stark contrast to the behaviour of general purpose, nonreversible integrators, for which the energy error increases linearly in time on a suitable time scale. Note that if the vector field f is Hamiltonian, then the integrator φ is close to symplectic. However, this plays no role in our analysis and we believe it is irrelevant. For any orbit {x(t) : t ∈ R}, we have that {R(x(−t))} is also an orbit. If it is the same orbit, it is said to be symmetric. Otherwise, it is said to be nonsymmetric. Any orbit that intersects the fixed set {x : R(x) = x} of R is symmetric. Symmetric orbits display typical conservative behaviour, for example, the eigenvalues of symmetric fixed points have the same symmetry as those of Hamiltonian systems, and there is a KAM theorem for symmetric quasiperiodic orbits [9, 15]. Nonsymmetric orbits, on the other hand, cannot ‘see’ the reversing symmetry and can have any dynamics, including asymptotically stable fixed points and strange attractors (whose image under R must be a strange repellor). This is commonly observed only in low-dimensional systems [9, 15, 6]. For typical high-dimensional systems such as those of molecular dynamics, the phase space consists of an ergodic ‘sea’ containing tiny islands of regular (e.g. quasiperiodic) orbits. Since there is no known mechanism which could keep a chaotic orbit in the sea bounded away from the fixed set of R, these orbits are believed to be generally symmetric. (If the system has nonsymmetric first integrals I : M → R , then only orbits starting on the fixed set {x: I (x) = I (R(x))} of I can possibly be symmetric. Therefore, we assume that the system either has no nonsymmetric integrals, or the system and integrator are both restricted to the fixed set of I. Momentum and angular momentum are examples of such nonsymmetric integrals.) By backward error analysis [5, 14], the integrator φτ is (exponentially close to) the time-τ flow of the modified vector field ẋ = f̃ (x) = f (x) + τfp(x) + O(τ) where f̃ and fp are R-reversible. Here p is the order of the method. Therefore, the energy evolves according to Ḣ =: h = if̃ dH = if dH + τpifp dH + O(τ) = τpifp dH + O(τ). H is symmetric by assumption, so Ḣ is antisymmetric: Ḣ ◦ R = Rif̃ dH = iR∗f̃ R∗ dH = i−f̃ d(R∗H) = −if̃ dH = −Ḣ . (1) Under these circumstances the evolution of H for a symmetric ergodic orbit can be modelled as a diffusion process as in [11]. The approximation is valid for time scales which are long enough that one can average over the fast motion in x but short enough that the total energy drift is small. On such an intermediate time scale we think of the orbit as consisting of a fast motion on an energy surface (which plays the role of the ‘angles’ for a diffusion process) and a slow drift in H (which plays the role of the action) transverse to this surface. Letter to the Editor L595 We now make the following key assumptions: (i) the flow of f̃ is ergodic on an invariant set A ⊂ M with respect to an invariant measure μ and (ii) the set A is symmetric. The measure is necessarily symmetric. Then, we have (recalling Ḣ = h) lim T→∞ H(T ) T = lim T→∞ 1 T ∫ T